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We investigate nonlinear stability of two classes of cosmological solutions in massive gravity: 
isotropic Friedmann-Lemaitre-Robertson- Walker (FLRW) solutions and anisotropic FLRW solu- 
tions. For this purpose we construct the linear cosmological perturbation theory around axisym- 
metric Bianchi type-I backgrounds. We then expand the background around the two classes of 
solutions, which are fixed points of the background evolution equation, and analyze linear perturba- 
^ ■ tions on top of it. This provides a consistent truncation of nonlinear perturbations around these fixed 

' point solutions and allows us to analyze nonlinear stability in a simple way. In particular, it is shown 

that isotropic FLRW solutions exhibit nonlinear ghost instability. On the other hand, anisotropic 
FLRW solutions are shown to be ghost-free for a range of parameters and initial conditions. 



I. INTRODUCTION 

O | Our current understanding of gravitation is based on general relativity (GR), whose predictions are consistent 
• with all available experimental and observational data. To this end, in many cases the so called parameterized 
' post-Newtonian (PPN) formalism has been adopted to quantify possible deviation from GR and to compare with 
experiments and observations. There are indeed many theoretical models that fall into the scope of the PPN formalism. 
' While the PPN formalism has been useful, however, it does not include finite range modification of GR. From 
phcnomcnological viewpoints, in analogy with the finite range nature of the weak interaction in the standard model of 
■"sj" ■ particle physics, one might argue that the possibility of finite range modification of GR should be taken into account 
| when experimental and observational data are used to constrain deviation from GR. Additionally, such an extension 
is also expected to alter behavior of gravity at scales larger than the Compton wavelength of the graviton; a graviton 
mass as small as the current expansion rate of the universe may provide the source for the accelerated expansion 
CO . alternative to dark energy. 

From the theoretical point of view, construction of a finite-range gravity theory has been a long-standing problem. 
The massive extension of GR at the linearized level was first considered by Fierz and Pauli |lj in 1939. It was then 
pointed out in 1970 that the linear theory suffers from a discontinuity @, Q in the massless limit, now known as the 
vDVZ discontinuity. The discontinuity can be alleviated by allowing for nonlinear terms a la Vainshtein Q but, as 
pointed out in 1972 by Boulware and Deser [f|, generic nonlinear massive theories contain an unwanted sixth degree, 
in addition to the five polarizations of massive spin-2 field, which leads to ghost instability. Since then, it took a 
long time until a nonlinear extension of the Fierz-Pauli theory without the Boulware-Deser ghost, dubbed the dRGT 
theory, was found recently 0, Q , where the unwanted extra degree is successfully eliminated by construction 0] . 
However, the absence of the sixth mode does not necessarily imply that the theory is healthy. For instance, in some 
backgrounds, one (or more) of the five polarizations of massive graviton may be superluminal jol— Tlll| . or may become 
a ghost fT2| . 

The analysis of the degrees of freedom in the dRGT theory, which originally uses the Minkowski fiducial metric, 
are done by considering the Minkowski background (see e.g.[7[), which resides in the "normal branch". On the other 
hand, for the theory to have a phcnomcnological relevance, it is essential to obtain non-trivial background solutions 
which can describe the cosmological evolution. However, in the normal branch no cosmological solution exists flsT ). 
If the theory is extended to allow for a general fiducial metric 1 , expanding universe solutions can be realized in the 
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1 The dRGT theory extended in this way was shown to be free of Boulware-Deser instability in Ref . [T3 | . 
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normal branch [15l - ll7| |. but the Higuchi bound [12( cannot be respected; for instance, for de Sitter reference metric, 
the scalar graviton becomes a ghost at expansion rates higher than the graviton mass fltl]. 2 

In the presence of a negative spatial curvature [2l| or a general reference metric (HI, there exists a second branch 
of solutions, dubbed the "self-accelerating branch" , which stems from the factorized form of the temporal Stiickelberg 
equation of motion (or equivalently, from the Bianchi identities). With the strict Friedmann-Lemaitre- Robertson- 
Walker (FLRW) symmetry, this branch is disconnected from the normal branch and as a result, the dynamics of 
perturbations are dramatically different; even though the massive theories should have five degrees of freedom after the 
removal of the Boulware-Deser degree, only two degrees (the two tensor modes in the standard SO(3) decomposition) 
propagate at the linear level . The resulting linear theory is thus very similar to GR, with the exception of a non- 
vanishing and time-dependent mass for the two gravitational wave polarizations, which may have left a characteristic 
signature on observables [22j ■ This peculiar behavior may be considered as a side-effect of the FLRW symmetry, which 
leads to the disconnection of the two branches. In fact, as we show in Sec. IIII1 the introduction of anisotropics to 
the background and thus deviation from the FLRW symmetry break the factorized form of the temporal Stiickelberg 
equation, leading to a coupling between the two branches. Around such deformed backgrounds, all five graviton 
polarizations are generically dynamical at linear order. 

From the technical point of view, one of the goals of the present paper is to construct the theory of cosmological 
perturbations around an axisymmetric Bianchi type-I universe in dRGT theory. From physical viewpoints, on the 
other hand, the goal of the paper is to investigate nonlinear stability of two distinct classes of cosmological solutions: 
self-accelerating isotropic FLRW solutions [HI, [2l[ and anisotropic FLRW solutions [23| . We will achieve this goal 
by allowing for small deviation of the background from the above mentioned two classes of solutions, applying the 
general formalism to analyze perturbations around those deformed backgrounds and then considering such linear 
perturbations as leading nonlinear perturbations around the undeformed backgrounds. 3 

The resulting quadratic kinetic action, which is non-zero for all five degrees at linear order in background deforma- 
tions, corresponds to the cubic kinetic terms from the perspective of undeformed backgrounds. 

This approach, previously reported in [25| for isotropic FLRW solutions, reveals that the isotropic FLRW solution 
always has a ghost mode at nonlinear order with a gapless dispersion relation. This conclusion is valid for all 
homogeneous and isotropic solutions in the self-accelerating branch. 

Although the deformation of the FLRW background with small anisotropy leads to an inevitable ghost mode, 
there may still be healthy regions with relatively large anisotropy. This is in analogy with scalar field models of 
ghost condensation (2(| [27| ■ As the second application of the formalism developed in the first half of the present 
paper, we thus explore this possibility, by considering the anisotropic FLRW attractor solution of Ref.[23| as the 
undeformed background solution. On this attractor, the expansion of the physical metric is isotropic, with a de Sitter 
evolution. The fiducial metric also exhibits a de Sitter expansion but is actually anisotropic from the viewpoint of 
comoving observers residing on the physical metric. This background thus contains anisotropy which can only be 
probed by perturbations, whose evolution depends on the alignment of the privileged direction. Applying the general 
perturbation theory to this example yields that on the attractor, there are still two degrees of freedom which do not 
propagate at linear order. By adding homogeneous deformations to the background, we find that the perturbations 
can in principle have a stable behavior at nonlinear order, depending on the parameters and initial conditions. This 
result suggests the existence of a healthy region between the isotropic and anisotropic FLRW solutions. 

The paper is organized as follows: we start by reviewing the massive gravity action in Scc|nJ Without choosing the 
solution, we first present the equations for an axisymmetric Bianchi-I background in Sec. IIII1 and metric perturbations 
in Sec lIVl In Sec|Vl we expand the background with respect to small anisotropy and discuss the implications regarding 
the nonlinear stability of the isotropic FLRW solution (this is a detailed calculation of the results already presented in 
[H||). In Scc lVIl we consider small deviations from the anisotropic fixed point solution of [23| and obtain conditions 
for linear and nonlinear stability of these solutions. Finally, in Sec lVIH we conclude the study with a discussion of 
our results. The paper is supplemented with several Appendices, where technical details are presented. 



II. ACTION 



We start this section by reviewing the action for nonlinear massive gravity in generic setup. 



2 An exception to this conclusion may be provided if gravity is partially massless, i.e. if there exists an additional symmetry which forces 
the Higuchi bound to be saturated, leading to the disappearance of the scalar mode within five propagating degrees of freedom f 1SH - 
However, there are indications that such a construction cannot be obtained from the dRGT theory llSL |20H - 

3 This is qualitatively similar to the instability encountered in the self-accelerating solutions at the decoupling limit, arising once the 
vectors are turned on [23 . 
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The graviton mass is introduced through a non-derivative term by [f| , 

I = I eh, a [g] + matter [g, Vj + ^massls" 1 /], (1) 
where the first term is the Einstein-Hilbert action with a cosmological constant A 

j\,t r 

Irh, a [g] = -P- / rf 4 ^f—g (R - 2 A) , (2) 



and the matter sector consists of fields {ip}, coupled to gravity through the physical metric g. The existence of 
a massive graviton implies breaking of the general coordinate invariance, which cannot be accomplished by the 
physical metric g alone Q. For this task, the space-time tensor f^, dubbed the "fiducial metric" is introduced and 
parametrized as 

f^ = I A B^ C )d^ A d^ B , (3) 

where the four scalar fields <f) A [A = 0, 1, 2, 3) are the Stuckelberg fields, which are responsible for the breaking of the 
four gauge degrees of freedom, while fAB(<t> C ) is the metric in the field space. 

Imposing the absence of the BD ghost in the decoupling limit, the mass part in |T]) can be constructed as [f| 

-;4- 



imass[<r /] = A£ m 2 I d 4 x yf=g~ (£ 2 + a 3 C 3 + a A U) , (4) 

with 



£2 = - 2 {n-[^]) , 



£3 = 6 



[/C] 3 -3[/C] [K 2 ] +2[/C 3 ]) 



U = y a ([icf -6[IC} 2 [K 2 ] + 3 [K 2 ] 2 + 8 [K] [/C 3 ] - 6 [/C 4 ] ) , (5) 



where the square brackets denote the trace operation and 



K$ = 6S- (y/F^f) ■ (6) 



III. ANISOTROPIC BACKGROUND 



In this section, we obtain the background field equations for an anisotropic extension of the FLRW physical metric. 4 
Since one of our goals with this extension is to obtain the stability conditions of the isotropic metric, in the isotropic 
limit, the whole system needs to reduce to the cosmological solutions in [TUHl]]. For this reason, we consider a fiducial 
metric that is of the flat FLRW form, 

ffiv = -n 2 ^ )^ ^ + a 2 (0°) (d^d^ 1 + Sijd^d^) , (7) 

where nonzero expectation value of the Stuckelberg fields coincide with the space-time coordinates ((f) ) = t, (tp 1 ) = 
x, (</>') = y t . In the rest of the paper, Greek indices span the space-time coordinates, while lower case latin indices 
i,j = 2,3 correspond to the coordinates on the y-z plane, with y 2 = y, y 3 = z. 

The physical metric is chosen to be the simplest anisotropic extension of FLRW, namely, the axisymmetric Bianchi 
type-I metric 

gf v dx* 1 dx v = -N 2 (t)dt 2 + a 2 {t) (e 4CT(t) dx 2 + e~ 2ff(t) 8 rj dy l dy 1 ^ . (8) 

Varying the Stuckelberg fields around the background values 

A = x A + n A , (9) 



4 For a more general discussion of Bianchi class A cosmologies in the context of bigravity, see [531 . 
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the mass term action up to first order in n a becomes 

/mass = 4°L + M|, m 2 J d 4 x N a 3 n n° £4, + O [(^ a ) 2 ] , 
giving the background equation of motion for the Stiickelberg fields as 

£$ = (H + 2E-H f e~ 2a X) + 2 J^ v) (H - S - H f e a X) = . 



where 



and 



H = 



aN 



H f = 



a n 



TV 



a 



(10) 



(11) 



(12) 



rO) - 



Ji y) = 3-e 



3- 2e CT X + a 3 (3-4e CT X 



a 4 (1 - e^AT) , 



la 



X - e a X + a 3 (3-2e~ 2 ' 7 X ~2e a X 



-a x 2 



a 4 (1 - e~ Ia X) (1 - e a X) 



(13) 



The expansion rate for the fiducial metric Hf is related to the invariants of the field space metric fabi'P ) which are 
fixed by the theory, and is independent of the choice of the background values of <j) a . Thus, Eq. fTTj) can be interpreted 
as an algebraic equation for a (or equivalently for X), but not a differential equation. We also note that, in the 

isotropic limit (a, £ — > 0), we have = , and the corresponding equation in the FLRW case (c.f. Eq.(A.ll) of 
Ref. [ll|) is recovered. 

Finally, we calculate the equations of motion for the metric g^. Since we are interested in the stability of the gravity 
sector only, we neglect any matter field and consider only the vacuum configuration. The independent components of 
the field equations for the physical metric are 



3 (H 2 - E 2 ) 



A 



3 - + 32TE 



6 f + 3S2 



- (6 + 4 a 3 + a 4 ) + (3 + 3 a 3 + a 4 ) (2 e a + e~ 2 a ) X 

-(1 + 2 a 3 + a 4 ) (e 2 a + 2 e~ a ) X 2 + (a 3 + a 4 ) X 3 

= m 2 (e- 2a - e a )X \(3 + 3a 3 + a 4 ) - (1 + 2a 3 + a 4 ) (e a + r) X + (a 3 + a 4 ) r e a X 2 



m 2 g X 



(3 + 3 a 3 + a 4 ) (2 e a + eT 2 a - 3 r) + (a 3 + a 4 ) 3 - (e 2 a + 2 e'*) r 



X' 



where we have introduced 



+2{l + 2a 3 + a 4 ) (2 e a + e~ 2 a ) r - (e 2 a + 2 e~ CT ) X 



NX 



, (14) 



(15) 



IV. PERTURBATIONS 



In this section, we calculate the action quadratic in perturbations around the metric ([8]). The most general set of 
perturbations around the axisymmetric Bianchi type-I are given by [29j 

( -2N 2 <$> ae 2a Nd x x ae~ a N (d,B + v,) \ 

g$=\ a 2 e 4 ^ a 2 e° d x (8,(3 + A,) , (16) 

\ a 2 er 2 ° [t 5 t] + d l 8 J E + d (l h j} ] J 

where duhj\ = (dihj + dihj)/2 and d l Vi = d l \ t = d l hi = 0. Note that, since the y-z plane is Euclidean, the indices 
i, j are raised and lowered with <5 y and <5y . Similarly, we decompose the perturbations of the Stiickelberg fields © as 



7T A = (7T°, dm\ aV + 7T l ) , 



(17) 
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where diir 1 = 0. Since the vectors are defined on the 2d y-z plane, the transversity condition can be used to reduce 
each of these vectors to a single degree of freedom 

v i = e i d j v i A, >/<l,\- h i = e i J d j h, tt; = e/dj 7r dd , (18) 

where is a unit anti-symmetric tensor with €23 = —£32 = 1 and = eikS ■* ■ 

Although the nonzero a in the background metric (|5J) breaks the isotropy, there is still a residual symmetry 
corresponding to rotations around the x axis. Under such a rotation, the perturbations which transform as 2d scalars, 
or even modes ($, Xi B, ipj P, T i E) and the perturbations which transform as 2d vectors, or odd modes (v, A, h ) 
decouple at the level of the quadratic action. This allows us to study each sector separately in the following. 



A. Gauge invariant variables 



Before moving on to the discussion of the action, let us first determine the transformation properties of the pertur- 
bations and build gauge invariant variables. We consider infinitesimal coordinate transformations 



with 



xfi _^ x m + ^ f 



Under this transformation, the even perturbations transform as 



<I> 

X 
D 

p 

T 

E 



while for the odd perturbations, we have 



Ne 



-la 



„2<7 



x - 

D 



a 

Ne °-e 



N 



a N 
if) + 2 N (H + 2 S) C° + 2 d 2 ^ 1 

T + 2JV(ff-S)(°, 



(19) 
(20) 



(21) 



v 

A 
h 

TTodd 



ae ■ 

v + —fj— ?odd , 

A + e- 3CT £ odd , 
h + 2 ^ oc jd , 

TTodd + £odd ■ 



(22) 



Using these transformations, we first define gauge invariant variables which do not refer to the Stiickclbcrg fields. For 
the even perturbations, we introduce 

1 



$ = $ d t 

2N 



H - £ 



X = X 
B = B 



2a(H -£) 
e CT 



ae a 

T — Of 



N 
a e~ 



-3(7 



-3<T 



2N 



E, 



2a(H- E) 

4> = ^-^§r-e~^d 2 x (2P^e-^E), 



(23) 
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while for the odd perturbations, we have 

ae~ CT ■ 

A = A- < ~y~ h ■ (24) 

The variables defined above are not sufficient to account for all of the physical degrees of freedom. Using the 
Stiickelberg perturbations, we can define three even variables 

fir = 7T° 



2N(H-T,) 1 



and one odd variable 



K = 7T dd ( 26 ) 

The gauge invariant variables defined in E qsJ23l) and (|24j) originate only from the perturbations of the physical metric, 
thus are already present in the GR case [301. The additional degrees in Eqs. ([25|) and (f26| are the physical degrees 
of freedom associated with the additional degrees of the massive graviton and arise from the breaking of general 
coordinate invariance by the nonzero expectation value of the Stiickelberg fields. 



B. Odd perturbations 



We start by the action for the odd perturbations only. The perturbed metric we consider is, from jSJ and (|T6 



ds: 



(id 



-N 2 dt 2 + 2 a e~ a N v, dt dy l + a 2 e ia dx 2 +2a 2 e a d x X l dx dy l + a 2 e 



2a 



5ij + d^jhj)) dy'dy 3 , 



while for the Stiickelberg fields, we consider 

l0 



= V 



(27) 



(28) 



After using these decompositions in the action dU, then switching to the gauge invariant variables defined in Eqs. (|24|) 
and (f26|) , the resulting action depends on the three perturbations (w, A, h^). Among these, v does not have any time 
derivatives and can be removed by solving the constraint equation. In general relativity, this operation also removes 
hi, and the final action can be written in terms of A only. In the nonlinear theory of massive gravity however, we 
expect that remains in the action as an extra degree of freedom coming from the Stiickelberg sector. 
We Fourier transform perturbations as 



S(t, x, y l ) 



dk r dk 



(2tt) 3 / 



T e i(k L x+ki y 1 ) 



5(t, k L , ki) , 



(29) 



where and ki are the components of the comoving momentum in the x and y l directions, respectively. Due to 
the 2d rotational symmetry around the x axis, the physical results are expected to depend only on the longitudinal 
component k^ and the magnitude of the transverse components kr = \/k\ + k 2 . 

Since the quadratic action does not depend on the derivatives of v, its equation of motion can be solved by 



ae" 



v = 



N 



p 2 {e a +r) +2e 2a m 2 X J ( 



(y) 



2m 2 g Xjf K 



■)|(e"A) 



(30) 



where is defined in Eq. (fT3|). while 



PL 



Pt 



2 i 2 

Pl+Pt 



(31) 
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are the physical momenta in longitudinal and transverse directions, and the total physical momentum, respectively. 
After a further field redefinition, 



Vi = -e 3<T A, V 2 



2e^pt 



X — 2h„. 



(32) 



the quadratic action takes the following form 

2 

.3 



iodd = ^Y J NdtdkL £kTaZ 



where 



and 



K n 



N 2 



K 



22 



7 Vi 



fiii IVil 2 - n| 2 |v 2 | 2 - n 2 2 (v?v 2 + v 2 *Vi) 



A n = 



K 



22 



2p 2 ' 

m 2 a 4 e- 2,J X J ( fp 2 T p 2 



4 2 m 2 g e 2 °X + (e CT + r) p 2 



n 2 n 



n 2 
"i 2 



F 1 

a 4 e _4,7 p|py 



to 2 e CT X 



2 2J 2 p|, 



1 + e 



-3(7 



[(1 + e- 3CT ) Ji -2J 2 ] 



m 2 a 4 Ap 2 p 4 
4p 2 (l + e 3CT ) 

m '° 4X ^ [2J 2P l + (l + e-^)JiP 2 T] , 



8(l + e 3ff ) 

Ji = (3 + 3 q 3 + «4) - (1 + 2 a 3 + a 4 ) (r + e" 2 ") I + (a 3 + a 4 ) r e" 2 " I 2 
J 2 = (3 + 3 a 3 + a 4 ) - (1 + 2 a 3 + a 4 ) (r + e CT ) X + (a 3 + a 4 ) r e a X 2 . 



(33) 



(34) 



(35) 



C. Even perturbations 

The perturbed metric for the even sector is, from ((8j> and (|16|) . 

rfs^vcn = -iV 2 (l + 2$) dt 2 + 2ae 2a Nd x xdtdx + 2ae~' T Nd l Bdtdy l 

+a 2 e ia (1 + V) cfe 2 + 2o 2 e ff 9 x 0i/3 efe dy i + a 2 e~ 2 a (1 + r) + d t d 3 E] dy l dy J , (36) 

while the even perturbations of Stiickclbcrg fields read 

cfP = t + tt° , (p 1 =x + d^ 1 , f =y l + 9V . (37) 

Applying these decompositions to the action (fl]), we then switch to the gauge invariant variables defined in Eqs. (f2"3")l 
and (|25[) . The resulting action is then manifestly gauge invariant and depends on seven perturbations ($, B, t/j, 

TV, /3,r, -E^-)- The three degrees arising from the <7g^ perturbations, namely S and x, are non- dynamical, thus can 
be integrated out. Furthermore, the kinetic term for the f w is proportional to background equations of motion, so 
this degree is also non-dynamical. We interpret this as the would-be BD ghost, which is eliminated in this theory by 
construction. In the massless theory (GR), using the constraint equations also removes the degrees (3 n , E^, leaving 
only ip in the action, which becomes one of the gravity wave polarizations in the isotropic limit. However, due to the 
nonzero mass of the graviton, these two degrees are generically dynamical. Thus, the metric perturbations in vacuum 
has three physical even perturbations. 
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Upon expanding the fields in momentum space, the resulting action is formally 
M 2 



/i 2 2„ = ^ / Ndtdk L d 2 k T a 3 



^K^~y^n 2 y + z^Ay + y f A T z + z^B^ + -^-B T z + z^cz 



, (38) 



with explicit expression for the matrices K, fl 2 , A, B and C given in Appendix [Bl In the above action, the basis 
y = (4>, $, En) T corresponds to the three dynamical degrees, while Z = ($, B, x, TVr) T contains the four non-dynamical 
ones. The equation of motion for the latter can be solved to give 



z = -c- 1 \^Ay + B^j , 

where we assumed detC ^ 5 . Using this solution in (|3"8"|) . we obtain 

M 2 



/i 2 L = ^ I Ndt dk L d 2 k T a 3 



_ v "vt - v 

— k — + — My + y^ m t — -y* n 2 y 

N N N N 



where 



K = K-B T C~ 1 B, M = —B T C~ 1 A . 



n 2 = n 2 + A T C~ 1 A. 



(39) 



(40) 



(41) 



V. NONLINEAR INSTABILITY OF FLRW SOLUTION 



In this section, we analyze the action for the perturbations studied in the previous section, by considering small 
anisotropic departure from a FLRW metric. As discussed in [251 ] . quadratic action in Bianchi-I background with small 
anisotropy, gives information on the nonlinear action in FLRW background. 

Assuming |oj <C 1 and \H/H\ <C 1, we introduce the parameter e <C 1 and expand the quadratic action with respect 
to e. 

The Stiickelberg equation of motion at leading order gives, 



implying 



3 + 3 a 3 + a 4 - 2 (1 + 2 a 3 + a 4 ) X + (a 3 + a A ) X 2 = 0(e) , 

jj* = -2 X [(1 + 2 a 3 + a±) — (a 3 + a 4 ) X] a + 0(e 2 ) , 

= X [(1 + 2 a 3 + an) - (aa + a A ) X] a + 0(e 2 ) , 
Ji = J 2 + 0(e) = (1 - r) X [(1 + 2 a 3 + a 4 ) - (a 3 + a 4 ) X] + 0(e) 



(42) 



(43) 



A. Odd perturbations 



Expanding the action quadratic in odd modes, the terms (|3"4"|) at leading order in the small anisotropy expansion 
arc obtained as 



K n = 



K 22 



4 ? 4 

2p 2 



0(0, 
<r + 0(e 2 ), 



4(l-r 2 ) 

7 = O(c), 
n 2 ! = K n x (p 2 + M 2 GW ) + 0(e) , 

Q 2 12 = 0(e), 

^ 2 = A- 22 x (c 2 p 2 )+0(e), 



(44) 



5 As we show in Sec lVIl this is a relevant assumption for the cases we study. 
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where 

1 - r 2 

M^ w = m 2 q {l-r)X 2 [(l + 2a3 + ai )-(a 3 +a4)X} , c 2 = — . (45) 

la 

Thus, at leading order, we identify the mode Vi with one of the gravity wave (GW) polarizations in FLRW background 
[lU . The extra degree of freedom V2 is massless and has sound speed c s . We also note that the leading order kinetic 
term for V2 may in general become negative, leading to a ghost degree. For |er| <C 1, in order to avoid the ghost and 
to make M GW positive, we require 

(l-r)a>0, (46) 
which also agrees with the condition for avoiding the gradient instability of V2, i.e. c 2 > 0. 



B. Even perturbations 



For the even sector, the components of the kinetic matrix defined in (|41|) at leading order are 



K12 
K 13 
K22 
K 23 



Pt 



4+ 0(e), 
■M 2 w pl(2p 2 +p 2 T ) 



2p 4 (l-r 2 ) 
a 2 M 2 w p 2 T (p 2 +pj) 



a + 0(e 2 ), 



4p 4 (l-r 2 ) 

9 n 4 M 2 n 2 
a a m GW P L 



a + 0(e 2 ), 



1 — r 2 

,4 71 r2 J2 J2 



a + 0(e 2 ). 



M 2 GW P 2 L P 2 T {p 4 r 2 -3H 2 M t 



GW, 



K. 



3H 2 p 4 (l-r 2 ) 2 
<^M GW p 2 T 



a 2 + 0(e 3 ) 



\-r- 



(47) 



Notice that at order e, the first mode has kinetic mixing with the second and third ones. One can rotate the basis to 
make the kinetic terms diagonal. The leading order terms of the eigenvalues of the matrix K can be found to be: 



K 2 



«3 



8p< 

2a i M 2 w p 2 L 
1-r 2 
a A M 2 w p 2 T 
1-r 2 



<r + 0(e 2 ). 



<7 + <D(e 2 



(48) 



The only eigenvalue which is non zero in the FLRW limit is the first one, and corresponds to one of the gravity wave 
polarizations. The remaining two are thus the extra degrees arising from the graviton mass. However, at leading 
order, the kinetic terms of these two modes are related through 



«3 



Pt 
"2p£ 



^2 , 



(49) 



which implies that for \a\ <ti 1, the two extra modes cannot simultaneously have positive kinetic terms; one of them 
is always a ghost. 

The dispersion relation for each mode can be obtained by a series of transformations, as described in detail in 
Appendix [Cl Assuming that the condition (|4l)|) is satisfied, the dispersion relations are found to be: 

lo\ = p 2 + M 2 w + 0(e), 

1 ^/ (10 p 2 + p 2 T f - 8 p 2 L p 2 T - {2 p 2 + 3 p 2 T ) 



oj. 7 = — 



24 a 
1-r 2 
24 a 



(10 p 2 +p 2 T f -8p 2 L p 2 T + (2p 2 + 3p 2 T ) 



(50) 
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The first mode corresponds to one of the GW polarizations in the FLRW limit, while the other two are extra degrees 
of freedom. Note that in this case, the second degree is a ghost. We also stress that dispersion relation for the extra 
degrees consists only of momentum terms and thus is gapless. 

VI. NONLINEAR STABILITY OF ANISOTROPIC FIXED POINT SOLUTION 

In this section, specializing the field space metric to be pure de Sitter (Hf = const.), we analyze the action for the 
perturbations studied in Sec llVl by considering small deviation from the anisotropic fixed point found in [23j 6 

H , <j 0} X constants, r = e~ 2o "° . (51) 

A. Expanding the background for small anisotropy 

Assuming \a — ao\ <C 1 and |£/i/o| <C 1, we expand the background quantities 

a = ao + ecn + 0(e 2 ) , 

£ = e£i+0(e 2 ), 

H = Ho + eHi+Oie 2 ), 

X = X + eX 1 +O(e 2 ), 

r = ro + en +C(e 2 ). (52) 

Note that Hf is not expanded since the fiducial metric is set to be de Sitter. Even if the fiducial metric were of general 
FLRW type, Hf would not be expanded as long as the unitary gauge is employed at the background level. 

Let us study the background equations order by order. At order 0(e ), upon using the characteristic relation for 
the anisotropic fixed point, ro = e~' 2<T °, or, 

X = |V-, (53) 

there are two more equations at this order; 

H 2 e 3a "{a 3 + a 4 ) - H H f (l + e 3ff0 )(l + 2a 3 + a 4 ) + H 2 f {3 + 3a 3 + a 4 ) = , 



3 + 77f e °(l + 2a 3 + a4 ) 



Hn to 2 

(1 + e 3ao )(3 + 3a 3 + a A ) +m 2 (6 + 4a 3 + a 4 ) + A = 0, (54) 



H 



f 



the solutions of which, determine the sets of fixed point values for H and ctq. However, noting that a 3 and a 4 enter 
these equations linearly, we can obtain a single set of solutions. Furthermore, introducing two mass scales as 



M 2 = -^-f [ff 2 e 3ffo (l + 2e 3CT »)(a3 + a4)-2Fotf/(l + e 3CTo + e 6CT °)(l + 2 a 3 + « 4 ) 
+F)(2 + e 3,T0 )(3 + 3a3 + a4)] , 



3 H TTl^ 

M 2 ee —3^ [iJ 2 e 6,To (a3 + a4)-2Ho^/e 3,To (l + 2a3 + a4) + ^(3 + 3a3 + a4)] , (55) 

which, along with the O(e') order background equations, allow us to express all physical background quantity in terms 
of Hq, (Tq, M and M, in mass units of Hf. Thus, from here on, we choose the Hf = 1 units. 



The global stability of this fixed point solution, shown in [23l |. can be interpreted as the validity of the cosmic no-hair conjecture [3l( in 
massive gravity. This point was addressed in the context of bigravity in 1331 . 
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Upon using the zero order relations above, the 0(e) equations yield, 
Hoe 2 " 



-Y, 



M 4 - 27 H 2 {3M 2 - M 2 ) 



2M 4 + 27 H%(3M 2 



2M Z ) 



a x + 18 Ho M Ei 



}■ 



M 



27H Q M 2 (7i +2M 2 Ei 



rie 2 ' = 



Si 



A/ 4 - 27 i?g (3 M 2 - M 2 ) 
243(3M 2 + 2M 2 )(3M 2 - M 2 )H$ + 9AI 2 H 2 (81M 4 - 66M 2 M 2 + 4M 4 ) - M 6 (27M 2 - 2M 2 ) 
[M 4 - 27H 2 (3M 2 - M 2 )][A~I 4 + 9H 2 (3M 2 - M 2 )} 
9H a (9M 2 - 2M 2 ) 
~il7 4 - 27i? 2 (3A/ 2 -M 2 ) 

3 H 



CTl 



M 2 M 2 (9# 2 + M 2 ) 



Af 4 + 9i/^(3M 2 -M 2 ) 



(Tl - Ei 



(56) 



Using these, we can express all order 0(e) physical quantity in terms of o\ and Ei only. 

Finally, we present the expansion of the J functions used extensively throughout the text. We use the above 
expansion procedure in Eqs.([T3|) and (|35|) . to obtain 



J. 



(■>■) 



2H 2 M 2 e 3a °(l 



O 3o- \2 



2(1- H ) 2 M 2 e 3CT o + 9 M 2 (l - H e 3 °o) 



0(e) 



Ay) 



■h = 



9H 2 M 2 e 3a °(l - e 3a °) 


3 M 2 (9 H 2 + M 2 )cti +2 H (9 AI 2 - 2 M 2 )E X 


M 4 — 27 .ff 2 (3 M 2 — M 2 ) 




2(1-H Q ) 2 M 2 e 3r7 ° +9M 2 (1- H e 3 °«) 2 



ou 2 



9H$M 2 (1 - e 3a °) 2 



0(e) 



2 (1 - H ) 2 M 2 e 3CT o + 9 M 2 (l - H e 3 °») 2 

27 HlM 2 M 2 e 3aa {\ - e 3CT °)(9 H o + M 2 )<n 



M 4 +9H 2 {3M 2 - M 2 ) 2(1- H Q ) 2 M 2 e 3<J ° + 9 M 2 (l - H Q e 3CT °) 2 



0(e 2 ). 



(57) 



B. Odd perturbations 



We now consider the odd sector action in Eqs. (|53")) - ([3"4"|) and expand the background quantities around the anisotropic 
fixed point solution, as prescribed in the previous subsection. The coefficients in the action now read, 



A n = 



fi 2 

"12 



2p 2 



0(e) 



K 22 = - 



3M 2 e 2a ° a 4 p 2 T 


3M 2 {9H 2 - 


V M 2 )a x +2H (9M 2 


- 2M 2 )Ei 


4(1- 


- e 6cr °) 


M 4 


- 27 H 2 (3 M 2 -M 2 ) 





12p|^E 



Pt 



nfi = Kn x [ p 



-+0(e 2 ), 
2 3A/ 2 p 2 



0(e), 



3Al 2 a 4 p 2 L p 4 T 

4 e 4 a ° p 2 
3 M 2 a 4 , p\* 



+ 0(e), 
0(e). 



0(e 2 ), 



(58) 



where the components of the physical momentum are now defined at order 0(e°), i.e. 
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After the small anisotropy expansion, we see that the first derivative mixing term 7AT22 ~ 0(e 2 ) and can be 
neglected. We introduce a new field basis 



Vi = y/K n Vi , 



V 2 



'K 2 2 V 2 ; 



(60) 



and bring the kinetic terms to canonical form. Note that we assume that there are no ghost degrees, i.e. Ku, K22 > 0. 
Since the contribution from the time derivatives of the rescaling to the mass term is suppressed in the e expansion, 
the odd sector action at leading order becomes 



I^i-^r / Ndtdk L d 2 k T a 3 



Vi 



N 2 



+ 



V, 



N 2 



n 2 



Vi - 



"12 



_(Kv 2 + v 2 *Vi)--^ |v 2 | 

/ K 11 y / K 2 2 -^22 



(61) 



It is then straightforward to find the leading order terms in the dispersion relations, 

P 2 



,1 



'0(e), 



M 2 p^(e- 6ao - 1) 


M 4 - 27 if, 2 (3 M 2 


- M 2 ) 




2M 2 


3 M 2 (9 H 2 + M 2 )a 1 +2 H (9 M 2 - 


-2M 2 )S X 



+ O(e ). 



(62) 



We see that the dispersion relations for small deviation from the fixed point, arc dominated by momentum terms. 
Therefore, although the kinetic term of one of the odd modes tends to vanish on the fixed point, we cannot, in general, 
integrate out this mode from the action since no mass- gap is present. 



C. Even perturbations 



The kinetic matrix K defined in (|41[) is considerably bulky for presentation. In terms of the leading order in small 
anisotropy expansion around the fixed point, the matrix is formally 



K 



I O(e ) O(e ) 0(e) \ 

O(e ) O(e ) 0(e) 
V 0(e) 0(e) 0(e) ) 



(63) 



Since -K"n -^22 — K 2 2 = 0(e°), there is only one mode with a vanishing kinetic term on the fixed point. The two 
0(e°) eigenvalues involve square- roots (since they are solutions to a quadratic equation) and are not suitable for 
presentation (nor calculation). However, going to a new basis with a non-orthogonal transformation, we can have 
(relatively) simpler eigenvalues. We introduce a transformation of the form: 



R 



Kl2 K33 — K13 K23 
Kl3 ^22—^12 K23 




1 

K33 







(64) 



^23— -^"22 ^"33 

which we use to define new basis y = Ry. In this basis, the kinetic matrix becomes: 



R 1 KR 



( 


«1 







( 







K 2 







\ 








K3 J 


\ 



( K n 



Kj 3 K 22 +Kf 2 K 33 -2 K 12 Kis K 23 
Kl 3 -K 22 K 33 

K 2 2 





K33 










^33 



(65) 



Notice that the eigenvalues Ki are not the same as the eigenvalues of the matrix K, since the transformation R is 
not orthogonal. However, R has unit determinant, so the determinant of the kinetic matrix does not change. This is 
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enough to study the signature of the eigenvalues. The eigenvalues can then be computed to be: 

n -1 

+ 0(e), 



H 2 



K3 



8M 4 



Pt M 4 + 9fi 2 (3M 2 -Af 2 ) 



2aj e 8a ° M 2 p\ 


9H 2 p 4 {M 2 


- 3 M 2 ) 


f M 4 p 2 L (-2p 2 +p 2 L ) 




M 2 p\(M 2 -3p 2 ) 2 - 9H 2 (M 2 - 


-3M 2 ) 


6p 4 + M 2 {-Ap 2 +p 2 L ) 



3M 2 e 2a ° a\p 2 T 


3 M 2 (9 Hq + M 2 )cti +2H (9 M 2 - 2 M 2 )£i 


(1- 


e 6<T °) 


M 4 - 27 H 2 (3 M 2 - M 2 ) 





0(e), 



0(e 2 



(66) 



The eigenvalues which are non-zero on the fixed point yield the following no-ghost condition on the parameters 

ki>0, k 2 >0. (67) 

The third eigenvalue, which is of order 0(e), is actually proportional to the kinetic term of a similar mode in the odd 
sector, with 



K3 = 4if 2 f + 0(e) 2 



(68) 



where -^22 d is given in the second line of Eq. ([55| . Satisfying the third no ghost condition K3 > is rather difficult, 
since it involves evolving quantities a\ and £1. 



D. Evading the ghosts 



Since we are looking for a region in the parameter space where none of the degrees of freedom are ghosts, we impose 
that the eigenvalues of the kinetic matrix in the even sector are positive, i.e. 



«i, «2, K3 > . 



(69) 



which ensures that the odd sector is also safe. 

For the first two kinetic terms in Eq. ([66[) , it is possible to find regions where a given mode with any momentum is 
not a ghost. The sufficient (but not necessary) conditions we adopt here are, 



~ 2 9 M 2 (9H 2 ~M 2 ) 
M 2 <0, M 2 < ^<0. 



(70) 



Thus, K\ and K2 can be made positive, regardless of the momentum of the modes. 
With a parameter set satisfying ([70]) . K3 > gives 



Sgn(l-e ff0 ) -3|AT| 9 H[ , - \M Z \ <n + 2 H 2 |AT | - 9 |AT | £1 



> 0. 



(71) 



which depends linearly on the deviation of the background from its fixed point, i.e. the two dynamical functions <j\ 
and Si. This means that regardless of the value of M 2 and M 2 , there may always be a region where <j\ and £1 
conspire to give a negative K3. On the other hand, close to the attractor, the evolution of these two quantities depend 
on each other. The equation of motion for o\ can be written as [23J 



£1 +3i? £i + M>i = 0, 



where 



Mi 



3 M 2 M 2 (9 H 2 + M 2 ) 



AI 4 + 9H 2 (3M 2 - M 2 ) 

First, the condition that the fixed point is stable, i.e. M 2 > 0, combined with the conditions ([70]) . yields 

9H 2 - I M 2 1 > 0. 



(72) 



(73) 



(74) 
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We also wish that a\ is over-damped, so that there are no out-of-control changes of signs: 

9H 2 >AM 2 . (75) 

The solutions to Eq. (fT2")) are easy to obtain, 

~-H t 



a = e 



(76) 



Eilt-,00 = ( "i^o + \]l Hi Ml . (77) 



At late times, the general solution approaches to the C2 = solution, where we can write 



3 „ 9 



So, in this regime, the condition (|71|) can in principle be satisfied by choosing the appropriate sign for o\. 
As an example, we choose the set of parameters used in [23| 



a 3 = -—, a 4 = l, A = 0, ^ = 20, (78) 



20 



implying 



As a result, we find, 



0.51, #0^ 15.35. (79) 



M 2 = -45079, M 2 = -460.66, V ^ ° - = -186.9, M 2 = 362.10, -i^ - = 168.10 , (80) 

27 Hq 4 



which satisfy the no- ghost conditions for the first two modes ([70]) , fixed point stability condition (|74|), and over- 
damping condition (J75J. For this example, the kinetic term for the third mode, along with the late time attractor 
assumption ([77)1 becomes 

k 3 ^ -41.85 <4z>| o-i . (81) 

Thus, the parameter set ([78]) does not lead to any ghost degree, provided that i. we are close to the late time attractor; 
ii. we are in the correct side of the fixed point, i.e. o\ < . 

VII. DISCUSSION 

Since the introduction of the dRGT theory, there has been several attempts to construct a cosmological solution 
that is fully consistent with the known expansion history of the universe. In the present paper, we have shown that 
the self-accelerating solutions introduced in [ll| [2l| have a ghost instability at nonlinear order. A technical source 
of this problem is that the kinetic terms of three among the five degrees of freedom are exactly proportional to the 
equation of motion of the temporal Stiickclbcrg field for the self-accelerating branch. These kinetic terms reappear 
in the cubic order, and their behavior can be understood by a deformation of the FLRW symmetries. These results 
were previously reported in (25j and are also compatible with the independent analysis of Ref. [33[ . 

The next step is to seek another class of cosmological solutions on which linear perturbations have non-zero (and 
positive) kinetic terms. To this end, one needs to ensure that the factorized form of the Stiickelberg equation of 
motion is broken, since kinetic terms of some perturbation degrees of freedom are typically proportional to one of the 
factors in this equation. Since the FLRW symmetry consists of homogeneity and isotropy, there are two possibilities: 
breaking cither homogeneity or isotropy. 

In the present paper we have considered the possibility of breaking the FLRW symmetry by introducing anisotropy. 
With a relatively large anisotropy, we found a healthy region with non-zero and positive kinetic terms for all five 
perturbation degrees of freedom. This healthy region is in a neighborhood of (but not exactly on) the de Sitter fixed 
point solution with an anisotropic hidden sector introduced previously in Ref. [23[. To be more precise, signatures of 
the resulting kinetic terms depend not only on the direction of deformation of the background from the fixed point, 
but also on the signature of its time derivative. If the initial condition is such that the evolution is close to the 
fixed point and that the deformation is in the correct direction, then all five graviton polarizations in the anisotropic 
background can have positive kinetic terms for a range of parameters. 
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Exactly on the fixed point, however, quadratic kinetic terms for two of the five degrees vanish, while all other three 
degrees have positive kinetic terms in a finite region in the parameter space. Hence, for a small deviation from the 
fixed point, even though appearance of the ghost degrees can be avoided, two degrees of freedom have very small 
(cubic order) kinetic terms and no mass gap (see Appendix [F]) . This is an indication that these degrees are strongly 
coupled and remains to be a source of concern. 

Nevertheless, the anisotropic FLRW solution studied here is the first calculable example of a stable cosmology in the 
dRGT theory of nonlinear massive gravity. One of technical advantages of this solution is that the spatial homogeneity 
and the SO (2) invariance of the axisymmetric background allows decoupling between even and odd sectors at the 
linear order. Moreover, the absence of ghost in a neighborhood of this solution may indicate the existence of region 
without ghost nor strong coupling somewhere between this fixed point solution and the isotropic FLRW solution. We 
also note that the perturbations around other Bianchi type universes remain unexplored. 

Alternatively, one may decide to break the FLRW symmetries by introducing inhomogeneities. For instance, there 
are such known solutions [H, I34l - l39j where the background dynamics is FLRW-like (or reduces to FLRW through 
cosmological Vainshtein mechanism [13|]) while the Stiickelberg sector is contaminated by spatial inhomogeneities. 
In principle, these may also evade the FLRW ghost, although in the lack of symmetry, the stability analysis of 
perturbations becomes technically challenging (although see [40[). 

A different approach to break the factorized form of the time Stiickelberg equation is to extend the theory, through 
introduction of new dynamics, e.g. by imposing a dilatation-like symmetry 14 II] or considering variations of the 
parameters with a scalar field |42| ■ Indeed, the resulting cosmological solutions [43l445| can still be fully homogeneous 
and isotropic, with non-zero kinetic terms for perturbations. The perturbation analysis of these extensions will be 
addressed in a separate work [46j | . 
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Appendix A: Nonexistence of less symmetric fixed point solution of Bianchi— I 

In this Appendix, we show that fixed point solutions of Bianchi-I type always have axisymmetry, as long as the 
expansion is constrained to be isotropic and the fiducial metric is de Sitter. We introduce the generic Bianchi-I metric 
as 



ds 2 



-N 2 dt 2 + a 2 



e 2(<r 1+ a 2 ) dx 2 + e -2 CT1 dy 2 + p -2a 2 ^ 



(Al) 



which, for o\ = cr 2 , a\ = —2 <J2 and 172 = —2 ci, reduces to the axisymmetric metric used in [23j . up to redefinition of 
coordinates. 

For the present discussion, the equations of motion for o\ and ct 2 are sufficient: 



1 

37V 
1 

3iV 



2 Si 



where 



771 

H (2 Si + S 2 ) = -f- ( e - ai - a2 



Si + 2S 2 +i?(Si +2S 2 ) 



JL ( e -"i—"i 



- e CT1 ) X [(3 + 3a 3 + a 4 ) - (1 + 2a 3 + a 4 ) (e a2 + r) X 

+ {a 3 + a 4 ) e" 72 r X 2 } , 

- e CT2 ) X [(3 + 3a 3 + a 4 ) - (1 + 2a 3 + a 4 ) {e ai + r) X 

+ (a 3 + an) e ai r X 2 ] , (A2) 



X = 



H = 



Si = 



S 2 = 



(T 2 



(A3) 



NX' "a 1 " aN' N ' N 

We now look for fixed points of the type found in |23l |. namely, which have isotropy in the expansion (Si = S 2 = 0), 
but not in the normalization (ci ^ 0, er 2 7^ 0). Under these conditions, Eq. (|A2j) can be written as 

( e -»i-oa _ e ai) [(3 + 3 a3 + a4 )-(i + 2Q' 3 + a4)(e <T2 +r)A + (a 3 + a4)e' T2 rA 2 ] = 0, (A4) 

( e -*i-*s _ e ^) [(3 + 3a 3 + a 4 ) - (1 + 2a 3 + a 4 ) (e ai + r) X + {a 3 + a 4 ) e ai r X 2 ] = 0. (A5) 
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Let us discuss the solutions to these equations. For Eq. (|A4[) , one of the solutions is 

e CT2 =e- 2CTl , (A6) 
implying axial symmetry around z direction. Similarly, Eq. (|A5[) has a solution 

e <7 I=e -2a 2) ( A7 ) 

which implies axial symmetry around y direction. Both of these solutions correspond to the solutions already studied 
in [23j]. Removing these, (|A4j) and ()A5|) become 

(3 + 3a 3 + 04) - (1 + 2a 3 + on) {e* 2 + r) X + (a 3 + a 4 ) e CT2 r X 2 = 0, 

(3 + 3a 3 + a 4 ) - (1 + 2a 3 + a 4 ) (e ffl + r) X + (a 3 + a 4 ) e 5l rX 2 = , (A8) 

for which, the solutions o\ and (T2 are the same, implying axial symmetry around x direction. The only alternative 
solution is X = (2 + a 3 )/[r (1 + a 3 )], which is realized only for a 4 = 1 + a 3 + a 2 . Moreover, the remaining equations 
require a further tuning between a 3 , m g and A, so we drop this solution as well. 

Thus, we conclude that all the de Sitter fixed points in Bianchi-I with isotropic expansion are axisymmetric. 



Appendix B: Explicit action for even modes around axisymmetric Bianchi— I 



In this Appendix, we present the explicit expression for the action quadratic in even modes. After switching to 
momentum space, the action has the form 



iilL = I Ndtdk L d 2 k T a 3 



y ] y 

— K — 

N N 



y j 'n 2 y + z j 'Ay + y^A T z + z^B^ + ^-B T z + z j 'cz 



, (Bl) 



where the two field vectors are 



y 



( * \ 



B 

X 



(B2) 



while K and fl 2 are 3x3 real symmetric matrices, C is a 4 x 4 real symmetric matrix, A and B are 4x3 real matrices. 
Up to background equations, the components of these matrices are 



K = m 2 g a 4 X 




















(B3) 



m 2 „a 4 e 2 ' J Xp 2 T r 



3 6 a 



{e 2 °r-l)J i 



(*) 



-,2 c 



(y) 






l+e" 



2 a 2 



e^p 2 L -e a p 2 L 



1+e- 3 ' \ 

2 >! / 



(B4) 



A 



( An A12 Au \ 

A21 



y A41 A42 A43 ) 



B 



/fti \ 

B21 B 23 

B 32 

V 641 B i2 B 43 ) 



^ Cn C12 Ci 3 Ci 4 \ 

C12 C22 C 23 C2 4 

C13 C 23 c 33 c 34 

y Ci 4 C2 4 c 3 4 C44 j 



(B5) 
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where 



Al2 

A21 

Ail 

A42 
A43 



2 2 2 a v t( x ) 2 

m a e XJ\'p L 



B 
2 2 

i g a , 



3m 2 NrXjji x) 
2 g (1 -e3^ ( H - H f*° X h 

3m 2 g a 2 e 4a N r X p 2 L ' 

1 - e 3fT 
m 2 a 2 e~° ' X rXp 2 T 



(H-H f e°X) 



J ( v ] (H-T,-H f rX) + J 1 (H - £ - H f e a X) + J 2 (H + 2 £ - H f e~ 2<T X)] (,B6) 



B11 
B21 

£32 
B41 

642 
£43 



F- £, 
,2 „3 . 



m~ a" e~" X Pt 



m 2 g a*e^XJ^p 2 L 

2 <T , 



1 + e l a r 



-m 2 g NrXJ { ; ] 



1 
2 

m z g a' e^NrXJ^'pl 

1 + e 2a r 
ml a 2 e~° ' N r X J { / ] p\ 



(B7) 
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Cu 
C12 

Cl3 
C14 

C22 

C23 

C24 

C33 

C34 
C44 



-6(# 2 -£ 2 ), 
ae- a p 2 T (2H + T,) 
2ae 2 >|(i/-E), 



-m 2 e- 2 CT TV r X 2 H f { j{ x) + 2 e 3 CT J, (y) 



1 2 2 / -2. 2 , 2 ^* J . 

1 2 (T 2 2 

m 2 ae~ a Nr 2 X v\ 



2 2 I 4(T 2 

■a p L [e p T 



2m g X J. 



(*)■ 



-2 cr 



+ r 



jM „2 t(w) 2 1 



-2(77-2) 



Ji + J, 



(y) 



r 

-2a 



(H + 2 E) ( 2 e CT J 2 + 2 e" 2 * 7 (J 2 - J^) - J 



H f e~ 2,J X 



2e 4a Ji + 2e CT J 2 



H f e- Za X{H + 2T l ) J), 



J. 



(x) 



2 J. 



(</) 



iV 



2 J 



(y) 



^ J N 



(B8) 



Appendix C: Diagonalizing the action for even modes near FLRW solutions 

In this Appendix, we analyze the quadratic action for even modes around a Bianchi -I background solution, em- 
ploying the small anisotropy expansion. We obtain a simple action for the even modes, through a series of field 
redefinitions and transformations. 



1. Applying small anisotropy expansion 

We start with the action (f^Oj) in the small anisotropy expansion as discussed in Sec IV Bl 



yl Kl A + K Ml y 1 + y\ Ml ^ - y\ n 2 y 1 

N N N N 



(CI) 



Since we will be applying several transformations, it is useful to keep track of each step, so we change the notation 
such that y — >• y± . K — > K\ , M —¥ Mi and fr — > f2 2 . The leading order terms of the matrix components are given 
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by 



(*i)n = + 

(^1)12 = 2 p 4 (l-r 2 ) CT + °( e ), 

f^l - * 2m gwPt(P 2 +Pl) a , gfc2l 



(^l)22 = _ 2a4 i M C T P| g + 0(£ 2 )i 



rRn ^ 4 ^P|P| ( P 4 r 2 -3ff 2 M 2 w .) 

(*Oaa = fl4 fj^ g + 0(£2) - (C2) 



W» = -^frS) (Vr + 34)<x + 0(e 2 ), 
(Ml)22 - ~ 3F(l-r 2 ) ff + 0(£ } ' 



3ff(l 

,2 J\/f2 „ 2 



(Ml)32 - 3ff(l-r 2 ^ fT + ° (e } ' 



Was - ^fz^- + ^ 2 )' ( C3 ) 
(n?) u = % + <^), 

{nj) 12 = l -a 2 M 2 GwP 2 L + 0{e), 

(^)i 3 = -\a 2 M 2 GwP 2 T + 0{e), 

(nl) 22 = \a±M 2 GW p\ (3p 2 + P 2 )+0( e ), 

(fi 2 ) 23 = ia 4 A4 H/ p 2 p 2 , + 0(e), 

(0 2 ) 33 = \ a 4 (3p 2 + j£) + 0(e) , (C4) 

Although the mixing matrix Mi can be made anti-symmetric by adding boundary terms, we keep it for now, since 
the following transformations will spoil its symmetries anyway. 
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2. Diagonalizing the kinetic matrix 



We start by diagonalizing the kinetic matrix. To achieve this, we apply the following field transformation, 



where the transformation matrix is 



y 2 = i?r x yi 



(C5) 



Ri 



1 

(Jfl)l2 (gl)33-(^l)l3 (Kl)23 

(K 1 yi 3 -(K 1 ) 22 (K 1 ) 33 

(-K~l)l3 (^l)22-(A:i)l2 (^l)23 
\ (K 1 )' A 2 3 -(K 1 )22(K 1 ) 33 (Kl) 3 3 / 



1 

(KlJ: 



(C6) 



and det(i?i) = 1. Wc stress that R\ is not orthogonal. The transformed action takes the form 



M 2 

l£L=-f I Ndtdk L d 2 k T a 3 



^^2^ + 4 m 2 y 2 + y\ Ml % - y\ n 2 y 2 



N 



N N 



N 



where the new matrices arc given by 

K 2 = Rj K\ Ri 
M 2 = RlM 1 R 1 +R(K 1 



E>T ^ Rl 

R l n r d dT ».rT R l R l „ R ^ 



N ' 



fi* = Rl ill R, - ^ Ah R, - Ri AI( -rr — -tt- K\ -rr 



N 



N N 



N 



(C7) 



(C8) 



For the vacuum configuration with A ^ 0, the order of time derives can be determined from the background equations 
in Sec lIIII as 



H = 0(e 2 ). 



PL 



-HNp L +0(e), p T = -HNp T + 0(e), 



A I, 



GW 



Using these, we can obtain the components of the new matrices as: 



Pt 



(*0ii = + 



(K 2 ) 22 

(^2)33 



2a i Al 2 w p 2 L 
l-r 2 
a 4 M 2 w p 2 T 



a + G(e 2 ), 



l-r 2 



a + 0(e 2 



-Mi 



GW 



2(1-0 



+ 0(e). (C9) 



(CIO) 



Pt 



8 Hp 2 



m gwP 2 lPt (2p 2 +Pt 



r 



{M 2 ) 21 = - 2 7- -r-o- + C>(e J 



2i/p 2 (1 - r) 
2 a 4 M 2 w p 4 r _ 
3 if (1 - r 2 

2 a iM GwPLPT r , , 

3if (1-r 2 ) 
4#p 2 (1 - r) 
3#(l-r 2 ) 

^ M GwPT r „ , /n/-^-, 

3i7(l-r 2 



( M2 ) 23 = - o ^f „ 2 y ^ + <^ 2 ), 



(^) 32 = o3T^r T+OCe 2 ), 



A/2 „4 



(ft 2 ) 12 = 0(e), 
(fi 3 ) 13 = 0(e), 



(fi 2 ) 22 = X -^Ml wV l (3p 2 +pi)+0(e), 
(0 2 ) 23 = ia 4 M^p|p 2 + 0(e), 
(ft 2 ) 33 = l^M^p 2 (3p 2 +p|) +0(e), 



3. Normalizing the kinetic matrix 

We now do the following field rescaling, 



y 3 = R 2 1 y2 



where 



R 2 = diag 



1 -v 2 



1 



2V2-V , 

pj, V 2cr ctMgwPl V a a 2 Mgw Pt 



l-r 2 



1 



where we assumed (1 — r)a > 0. With this transformation, the action becomes 
M? 



4vL = ~Y Ndtdk L d 2 k T a 3 



¥ Ks 17 + f A h y3 + y * M * ¥ _ ^ 



where the new matrices are 



K 3 = R 2 r K 2 R 2 , 



M 3 



i2£ M 2 i? 2 + ^ # 2 



^2 



Rl2 MR pT 71 rT ^2 -R 2 „ #2 

— M 2 R 2 -R 2 M 2 —- W K 2 —, 



22 



with components, 



( l + 0(e) 0(e 3 / 2 ) 0(e 3 / 2 ) \ 

0(e 3 / 2 ) -l + 0(e) 0(e 3 / 2 ) 
V 0(e 3 / 2 ) G(e 3 / 2 ) l + 0(e) J 



(C17) 



(M 3 ) n 

(M 3 ) ia = (M 3 ) 13 - (M 3 ) 21 = (M 3 ) 31 

(M 3 ) 22 

(M 3 ) 23 = - (M 3 ) 32 

(^3) 33 



El 

H 



0{e) 



= -H 



3H 
\/2p L p T r 



3H 
3H 



E 
2~ 



iV(l + r) 
0(e), 



E 
2^ 



0(e) ; 



iV(l + r) 



0(e) 



(C18) 



(n§) 
(«!) 

(n§) 



= M 



13 



22 



2:i 



33 



GW 



0(e), 



(l-r 2 )(3p 2 +p 2 L ) | g(c0) 



12 cr 

(1 - r 2 )p L p T 
6V2a 



O(e ). 



(l-r 2 )(3p 2 +p 2 T ) | p Q 

6(7 



(C19) 



4. Anti-symmetrizing the mixing matrix 



We now add the boundary term: 



a?yi (M 3 + Mi ) y 3 



giving the new action, 



Mi 



I&L = ^2. I Ndtdk L d 2 k T a 3 



yV ,. 3> 3 , yl 



K s^7 + ^7 M i 3^3 - yl Mi ^ - yl n 2 y 3 



where the mixing matrix is now antisymmetric 



M 4 = - (M 3 - M 3 T ) = 



N N N 
( 0(e) O(Vi) 



O(v^) \ 



f^+C(e) 0(e), / 



and the frequency matrix is 

»4-"3+ 2 JVo3 dt 



a 3 (M 3 + M 3 T )] 



/ p 2 + M 2 w + 0(e) 



(C20) 



(C21) 



(C22) 



O(V^) 
jp 2 +pj) . 



(l-r 2 )(3p 2 +p 2 ) | g ^Qj (l-r 2 )p L p T 



12 cr ' ^ v " y 6 V2 (T 

(l-r 2 ) Plj p T 1 /n ^Q\ (1-r 2 ) (3p 2 +p|Q 



0(e°) 



(C23) 
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At this point, we see that the first mode is decoupled at the leading order, and has the same dispersion relation as 
GW in FLRW, 

ojf= P 2 + M% w + 0{e). (C24) 
From here on, we concentrate on the remaining two degrees only. 



5. Removing the mixing 



The action for the two degrees is now, 
r(2) Ml 



even, 2, 3 



— ^ / Ndtdk L d 2 k T a 
where y 5 = [(y 3 ) 2 , (^3)3] and 



y^ y\ + 



-l + 0(e) 0(e 3 / 2 ) 
0(e 3 / 2 ) l + 0(e) 



(C25) 



(C26) 



M 5 = P ^Z | ° 1 \ +0(e) 



2V2H \ i 0, 



(C27) 



m = 



The next field transformation is 



where 



(1 -r 2 ) / (3p 2 +p 2 L ) V2 P lPt 
12 * \ V2PLPT 2(3p 2 +p 2 T ) 

y 6 = R5 1 y 5 , 



Rn = 



where the transformation function satisfies 



In the new basis, the mixing matrix becomes 



cosh[0(i)] sinh[0(i)] 
sinh[0(i)] cosh[0(i)] 

0_ PL Ptt 
N~ 2s/2H' 



M 6 = Rl M 5 R 5 + R% K 5 ^f- = O(Ve) ; 



N 



O(e ). 



whereas the frequency matrix reads 



.^M 5 R 5 -R^M^-^K 5 ^- = R^lR 5 + O(e ), 



N 



Finally, the action is, at the relevant order 

M 2 



r(2) 

even, 2, 3 2 



j Ndtdk L d 2 k T a 3 



y^K^-yURlnjRM 



(C28) 



(C29) 



(C30) 



(C31) 



(C32) 



(C33) 



(C34) 
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6. Eigenfrequencies 

At this point, to diagonalize the frequency matrix, we still need time dependent transformations, which will in- 
evitably reintroduce the mixing. In that sense, it does not seem possible to diagonalize the system at the level of the 
Lagrangian. This is very similar to the coupled bosons in external background, where the diagonalization is done at 
the Hamiltonian level [47| . The extended calculation to include ghost degrees is presented in Appendix Q21 where the 
eigenvalues of the frequency matrix turn out to be the same as the eigenfrequencies in the diagonalized Hamiltonian. 

Thus, the eigenvalues of Q§ will actually give us the correct dispersion relations. We diagonalize the matrix Qq 
through, 



where 



Rl n 6 R 6 = Rl Rl n 5 R 5 R 6 = diag i 



-OJn 



!). 



(C35) 



Ra 



cosh[ip(t)} smh[ip(t)} 
smh[ip(t)] cosh[ip(t)} 



(C36) 



and ip satisfies 

cosh [2 (6 + ip)} 



10 p 2 +p 2 T 



As a result, the eigenfrequencies are found to be 7 

,.2 ~ 



L0 = — 



LO 1 = 



X-r 1 
24 a 
1-r 2 
24 a 



sinh [2 (0 + ip)}=- 



2,/2 PL p T 



'(10 p 2 + P 2 T f~ 8 V \v\ 



(C37) 



(10 P 2+P 2 T ) 2 -8p 2 L p 2 T - (2p 2 + 3p 2 T ) 



(10p 2 + P 2 T f - 8 P 2 L p 2 T + (2p 2 + 3p 2 T ) 



(C39) 



Appendix D: Dispersion relations of coupled system with ghosts, in external background 



The action of the modes in Section IV Bl contains a time dependent and non-diagonal frequency matrix, and one 
of the degrees has a negative kinetic term. At the level of the Lagrangian, it is not possible to recover a diagonal 
form. In this Appendix, we extend the formalism of [4?} to include ghost degrees, show that the Hamiltonian can be 
diagonalized and obtain dispersion relations of each degree of freedom. 

We start by an action of the form 



/ 



d 3 kdtC k = I I d A kdt 



0+ (t, k) r, fa k) - 0t ( tj k ) n 2 (t, k) k) 



(Dl) 



where <f> is a N = Ni + N2 dimensional array of fields, with N\ ghosts and N2 physical fields. f2 2 is a N x N real and 
symmetric matrix, which stays invariant under k — > — k. The signature of the kinetic matrix is 



-1 
+1 



(D2) 



7 We remind that in this calculation, we assumed that the second mode is a ghost, i.e. (1 — r) a > 0. Assuming that the third mode is a 
ghost (1 — r) a < 0, one obtains similar expressions for the eigenfrequencies: 



U 2 



(^Pj V( 10 P 2 +Pt) 2 ~ 8 P 2 lPt ~ ( 2 P 2 + 3 Pt) 



24 a 
r 2 -] 
24 a 



'(lOp 2 +p 2 T ) 2 - 8p 2 L p| + (2p 2 + 3p%) 



(C38) 
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The conjugate momentum array can be immediately found as 



— 7T- = 0J , 71^ = — f - = <\>) T}ji 



giving the Hamiltonian 



H 



d 3 kH h 



1 



7T f (i , k) tt(* , k) + ^ (t, k) ft 2 (t, k) </>(/; , k) 



The Hamiltonian equations of motion are: 

= 77 7T , 7r = — n 2 . 
We also define the matrix C = C(t, k) which diagonalizes il 2 through 

C T n 2 C = r 1 uj 2 (diagonal), 
where C is an element of group SO(N2, Ni), i.e. 

C T rjC = r}. 



(D3) 

(D4) 

(D5) 
(D6) 
(D7) 



We also note that C is invariant under k — > k. 

It turns out to be useful to include the matrix C in the decomposition of the fields <f> an d momenta 7r into mode 
functions in a basis of N dimensional creation/annihilation operator arrays a and as 



<f>(t,k) = C(t,k) [/i(t,k)o k + ft*(t,k)o^ 



7r(/j,k) = [C* T 0,k)] 1 h{t,k)a k + h*(t,k)a[ 



where mode functions h and their conjugates h arc N x N matrices. With this decomposition, we can rewrite Eqs. (|D5 



(D8) 



h = rjh — Th, 



h 



-r]uj 2 h + T T h, 



(D9) 



where T = C 1 C corresponds to the rate of change of the transformation C. 
Using the decomposition (|D8|) . the Hamiltonian density becomes 



h f rjh + h f t]lj 2 h & rj h* + rj oj 2 h* 
~h T rj h + h T uj 2 h h T rjh* + h T ij oj 2 h* 




(D10) 



Finally, we do a further redefinition, based on the solutions to the equations of motion in adiabatic regime, 



h = 



<2lu 



(a + 0) , 



(Dll) 



where a and (3 are N x N complex matrices and are generalizations of Bogolyubov coefficients. We stress that this 
redefinition keeps generality. 

Using the definition (|D11[) in the Hamiltonian density (|D10|) . we obtain the Hamiltonian density 



n k = 



where the new creation/annihilation operators are defined as 




(D12) 




(D13) 
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Since "qui is already diagonal, the Hamiltonian in the new basis is also diagonal. After normal ordering, we end up 
with 

1 N 

»=i 

t=l i=l 

This calculation shows that the eigenvalues of the frequency matrix in Eq. (|Dl[) actually correspond to the eigen- 
frequencies of the modes in the diagonal basis. 

Appendix E: The situation of perturbations on the fixed point 

Since the even mode action in Sec lVI Cl is very bulky, it was not possible to obtain the dispersion relations for each 
degree of freedom, at least not for arbitrary momentum. Observing that the matrix M which mixes first derivatives 
with the fields is of order O(e ), it is then a fair question to ask whether the kinetic term is modified when the action 
is brought to a diagonal form. 

However, this reasoning requires integrating out the field 3^3 = E x , a procedure which cannot be performed in 
our case as the field 3^3 does not have any mass-gap. Let us consider this issue in more detail. The term p, 2 in the 
Lagrangian, where C 3 — n 2 y$, is surely of order O(e ) but multiplied by a factor plj,. Instead, the kinetic term of 
this same field 3^3 ~~ which tends to exactly vanish on the fixed point - is multiplied only by a factor pq, as shown 
in Eq. (|66[) . Rescaling the field by 3^3 — ► y^/pr, makes the C(e)-kinetic term momentum-independent; however, its 
mass - being proportional to p\ - still tends to vanish for small pr- This is tantamount to saying that we cannot in 
fact describe the mass of the field in the small transverse-momentum limit, that is pt — > 0, because we still have that 
lim PT _ i .o(M 2 /Pt) ~~ ^ and, consequently, we cannot integrate the mode 3^3 any more in this limit. 

Therefore, we are led to deduce that we cannot, in general, study the perturbations on the exact-fixed-point solution, 
as this would lead to the inconsistency of being able to integrate out a mode which actually cannot be integrated out 
- at least for small pr- On the other hand, it makes sense to study that background which is not the fixed point, but 
close enough to it. We have followed this last procedure in our analysis in this paper. 

This situation also implies that, in order to check whether there is a strong-coupling limit on the fixed point solution, 
we need to analyze the Lagrangian at higher orders (cubic, at least) in the perturbations, and check whether there is 
indeed any non-trivial strong contribution/back- reaction coming from these higher order terms. 

Appendix F: Even sector dispersion relations in the IR 

Unfortunately, for the perturbation analysis around the anisotropic fixed point discussed in Sec lVIl the problem 
is technically very involved, and calculation of the full dispersion relations is very difficult. However, it is possible 
to finalize the diagonalization in the IR regime. This approach allows us to determine whether the even mode with 
vanishing kinetic term on the fixed point suffers from strong coupling or not. Notice that we already have one mode 
(in the odd sector) which has vanishing kinetic term and has no mass gap. 

Due to the direction dependence of the background, in the IR regime where the momentum dependence is removed, 
the dispersion relations still depend on the orientation of the momentum vector. We decompose the momenta as 

PL=p£, Pt=P^1-£, 2 , (Fl) 

and take the limit p 0. Proceeding the same way as we did for the FLRW solution (see Appendix [Cl for details), 
the diagonalization in this limit reveals, 

wi = O(p ), u 2 =0(p°), w 3 Uo = 0> (F2) 

where modes labeled 1 and 2 correspond to degrees which have non vanishing kinetic terms on the anisotropic fixed 
point, while for mode 3, the kinetic term is zero. 

Combining our result from the 2d vector sector, we observe that: 

• Even though by appropriate choice of parameters, we may avoid ghost degrees, for the two modes which have 
vanishing kinetic terms, there is no mass gap. Therefore, we expect that they are infinitely strongly coupled. 
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• The mode with 0(1) kinetic term in the odd sector is massless on the fixed point, with sound speed c s = 1. 

• The two modes with 0(1) kinetic term in the even sector have masses, but these depend on the orientation of 
the momentum vector. The gradient part of the dispersion relation was not calculated. 

For the two massive modes in the even sector, the mass terms are still to bulky for presentation. We end this 
section by looking at two extreme examples for the momentum direction. 



1. £ — > 1, momentum is in the x direction 

In this case, the momentum is aligned with the privileged direction x. The masses for this situation is 

2 3 M 2 M 2 (9 H'q + M 2 ) 



-3Ar , 



9 H 2 {M 2 - 3 M 2 )-M 4 



(F3) 



If the no-ghost condition (|70|) and fixed point stability condition (|74[) are satisfied, both squared-masses are positive. 



2. £ — !> 0, momentum is perpendicular to the x direction 

In this case, the momentum is along the y — z plane. The masses of the modes then become 



M 2 M 2 



M 8 + 36 Hi M 4 (M 2 — 3 M 2 ) + 243 H$ (M 2 - 3 M 



2 \2 



M 2 M s + 27 Hq (M 2 - 3 M 2 ) 2 (3 M 2 + 2 M 2 ) + 18 Hi M 4 (3 M 4 — 4 M 2 M 2 + M 4 ) 



M 4 



9 H'q (M 2 — 3 M 2 ) - M 4 



54 Hi M 2 M 2 



4 Hi ( M 2 — 3 M 2 



9H 2 (M 2 -3M 2 ) - M 4 



1 - 



9Hl(M 2 -3M 2 ) - M 



12 Hl{M 2 -3 M 2 ) 

Again, once the conditions ([70)) and (|74|) are imposed, both squared-masses are positive. 
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